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Abstract. In this paper we analyze a suite of cosmological simulations of modified gravitational 
action f{R) models, where cosmic acceleration is induced by a scalar field that acts as a fifth force on 
1-H all forms of matter. In particular, we focus on the bispectrum of the dark matter density field on mildly 

non-linear scales. For models with the same initial power spectrum, the dark matter bispectrum shows 
significant differences for cases where the final dark matter power spectrum also differs. Given the 
0^ different dependence on bias of the galaxy power spectrum and bispectrum, bispectrum measurements 

can close the loophole of galaxy bias hiding differences in the power spectrum. Alternatively, changes 
^ in the initial power spectrum can also hide differences. By constructing ACDM models with very 

• • similar final non-linear power spectra, we show that the differences in the bispectrum are reduced 

^ ^ (< 4%) and are comparable with differences in the imperfectly matched power spectra. These results 

indicate that the bispectrum depends mainly on the power spectrum and less sensitively on the 
$H gravitational signatures of the f{R) model. This weak dependence of the matter bispectrum on 

^ gravity makes it useful for breaking degeneracies associated with galaxy bias, even for models beyond 

general relativity. — 
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1 Introduction 

Observations of Type la supernovae suggest that the Universe has been accelerating since redshift 
z ~ 0.5 [11, 12]. Today the physical mechanism responsible for this process is still a mystery. The 
simplest model to explain the acceleration of the Universe is the ACDM (Lambda Cold Dark Matter 
model). This model assumes that the acceleration is driven by an exotic form of energy with negative 
pressure that might be related to the vacuum energy of quantum field theories. This theory is 
equivalent to adding an integration constant to the Einstein equations. 

Alternative theories to the vacuum energy propose a modification of gravity in the infrared 
that would produce an accelerated expansion. One possibility are the f{R) class of models (see [19] 
and references therein). These models produce accelerated expansion through a modification of the 
Einstein-Hilbcrt action by an arbitrary function of the Ricci scalar R. As a consequence, an extra 
propagating scalar field appears that mediates a fifth force on all forms of matter. The range of this 
force depends on the functional form of f{R). In order to satisfy solar system tests, f{R) models are 
often chosen to present a chameleon behavior. The chameleon mechanism makes the extra scalar field 
become increasingly massive in higher-curvature regions, suppressing the range of the fifth force in 
dense environments. 

In previous works, cosmological simulations [9] have been used to study the power spectrum 
[10] and halo statistics [13] of these kinds of models. More recent studies with higher resolution 
have confirmed these previous results [21] and extended the investigation to smaller scales. In the 
present work we focus on how the dark matter bispectrum is modified in this class of models. While 
these models also predict a non-linear matter power spectrum different from the ACDM one, it is 
nevertheless interesting to look at the bispectrum for at least two reasons: a) except for gravitational 
lensing, measurements of clustering yield the galaxy or the baryon power spectrum, not the dark 
matter one: as baryonic physics and galaxy formation are complicated phenomena, the observed 
power spectrum may be biased, i.e. may differ significantly from the dark matter one; the bispectrum 
is well known for helping disentangle effects of gravity from effect of biasing e.g., [4, 20]. b) once we 
allow ourselves to consider non-standard models, the initial (linear) matter power spectrum does not 
have to be the power-law ACDM one to reproduce the observations. The form of the bispectrum 
kernel is a possible "signature" of gravity as it gets modified by any modifications from GR behavior 
e.g., [15]. 
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Here we pay special attention to see whether the bispectrum can be used to break degeneracies 
between models with the same observed power spectrum and the same cosmology, but different gravity. 
We begin in §2 with a review of non-linear gravitational dynamics in f{R) models, in §3 we briefly 
describe the simulations and in §4 we introduce the density field statistics. We discuss the results in 
§5 and conclude in §6. 



2 f{R) Gravity 

The f{R) class of models generalizes the Einstein-Hilbert action to include a function f{R) of the 
Ricci scalar R, 

'R + fiR) 



S = d X 



167rG 



(2.1) 



Here is the Lagrangian of matter and we have assumed c = h = 1. For standard GR with 
a cosmological constant, f{R) = —IGnGpA, whereas for modified gravity, the force modification is 
associated with an additional scalar degree of freedom = df /dR. In particular, in this paper we 
use the model for f{R) proposed by [6], 

where A is a constant with dimensions of length squared. We can write this equation as a function 
of its derivative evaluated at Rq (the background curvature today), namely fuo- We adjust the 
proportionality constant to match some effective cosmological constant p\ in the limit where fjio — ^ 0. 
For high enough curvature such that AR S> 1, f{R) can then be approximated as, 

/(i?) = -16^GpA-/flo^. (2.3) 

The modified Einstein equations can be computed by varying the Einstein-Hilbert action (Eq. 2.1) 
with respect to the metric. We work in the quasistatic limit where the time derivatives are negligible 
compared to the spatial derivatives. In this regime, valid on scales much smaller than the horizon 
1/H, the trace of the modified Einstein equations yields the fu field equation, 

y^Sfn = y [SRifii) - SnGSp^] , (2.4) 

where a is the scale factor, Sfn = fniR) — fniR), 5R ^ R — R and Spm — Pm — Pm- Here R is the 
background curvature that can be approximated by a ACDM universe for |/_r,o| ^ 1 ^nd pm (Pm) is 
the (background) matter density. 

On the other hand, the time-time component of the Einstein equations yields the modified Poisson 
equation, 

_o IGttG 9 „ „ , , , , 

V^vl/ = —a^Sp^ - -SRifn) (2.5) 

where ^ = SgoQ/{2gQo) is the Newtonian potential. 

For small fluctuations of the field, we can approximate 5R ~ {dR/dfii)\j^6fu. We will refer to 
this linearization as the non-chameleon limit. Conversely if the field fluctuations are large enough 
such that (5i?(/jj) cannot be linearized, the chameleon mechanism operates. We will refer to use of 
the exact, as opposed to linearized, equations as full f{R) or just chameleon models. 

The linearized field equations formed by Eqs. 2.4 - 2.5 can be solved for the Newtonian potential 
as a function of the density field. In the linear approximation for SR, these two equations in Fourier 
space yield, 

'4 1 fi^a^ 



k^^ik) = -A^G [ - - 3 ^,;_,^, ) a^5pM. (2.6) 
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This equation is identical to the one in GR but with a modification of the gravitational constant, 

4 1 fl{t)Mtf \ .9 7^ 

3 3 P + /i(i)2a(t)2 )■ ^ -'^ 

Here fj.{R) = {MfR/dR)-'^''^ is the effective mass of the scalar field fa and ^ just stands for n{R). The 
dependence on time is introduced though R{t). Note that when — 0, G^e G and we recover the 
ACDM limit, as expected. It is interesting to see that for a given value of fuo there are two different 
regimes for GeS, depending on whether the physical scale we are studying is larger or smaller than 
the inverse mass of the field. On large scales k ^ fi{t)a{t), GcS G and gravity behaves as GR, 
whereas on small scales k » ii{t)a{t) GeS 4G/3 and gravity is stronger than in GR by a factor of 
4/3. 

In other words, in Eq. 2.6, one assumes that the mass of the scalar field /i only depends on 
time and is the same in all regions of the Universe at a given epoch. However, for cosmologically 
interesting values of ^ the field is then essentially massless within the Solar System. The presence 
of such a scalar field (fifth force) is ruled out by light deflection and time delay measurements in 

— 1/2 

the Solar System, which are all consistent with GR. In the full non-linear f{R) theory, R (x ff^ 
can become very large in dense environments, suppressing the field and restoring the GR relation 
6R = SnGSpm (Eq. 2.4). Thus, gravity is not modified in the same way everywhere, but depends on 
environment. In regions with large potential wells (inside halos) the mass of the scalar field becomes 
large and therefore the effective range of interaction of this field shrinks recovering GR. We call this 
the chameleon mechanism. 



Geft(k,0 =G 



3 Simulations 

The simulations used in this paper are described in previous works [9, 10, 13]. Briefly, the field 
equation for fn (Eq. 2.4) is solved on a regular grid using relaxation techniques and multigrid iter- 
ation. The potential is computed from the density and ffj fields following Eq. 2.5 using the fast 
Fourier transform method. The dark matter particles are then moved according to the gradient of 
the computed potential, — V^*, using a second order accurate leap-frog integrator. 

The simulations were run using the values of Ifml — 10^*, 10"^, 10~^ (both for chameleon 
and non-chameleon cases) and 0, which is equivalent to ACDM^. The background expansion history 
for all cases differ from ACDM only at 0{fRo) and are hence practically indistinguishable. The 
cosmology used is Ha — 0.76, ftm = 0.24, = 0.04181, Hq = 73km/s/Mpc and initial power in 
curvature fluctuations = (4.89 x 10"^)^ at fc = 0.05 Mpc"^ with a tih of = 0.958. This initial 
power spectrum does not include the effects of baryon acoustic oscillations. Specifically, the initial 
conditions for the simulations were created using ENZO [8], a publicly available cosmological N-body 
+ hydrodynamics code. ENZO uses the Zel'dovich approximation to displace particles on a uniform 
grid according to the initial power spectrum. In order to propagate the initial power spectrum until 
late times the transfer function from Eisenstein & Hu [3] was used. 

The simulations were started at a = 0.02 and are integrated in time in steps of Aa — 0.002. 
All simulations used here correspond to boxes of comoving size L = 256 and 400 Mpc/h with 512'^ 
grid cells and 256'^ particles. For each box size, we have 6 runs for each value of fno, with different 
realizations of the initial conditions. 

4 Power Spectrum and Bispectrum 

The simplest statistic of interest of the matter density field is the power spectrum P(k), defined by 
the second moment of the Fourier amplitude of the density contrast, 

(5(k)<5(k')) ^ (27r)35^(k + k')F(fc) , (4.1) 

-"^Iii this paper wc arc always using negative values for /ijo, so when we talk about /^q we refer to its absolute value. 
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where (...) denotes the ensemble average over different realizations of the Universe. By statistical 
isotropy, the power spectrum does not depend on the direction of the k-vector. In practice we only 
have one observable Universe, so the average (...) cannot be computed. However, using the isotropy 
of the power spectrum we can compute the average over all different directions for each k-vector. 
Note also that P{k) is defined to be real. Since k = — k', (5(k)(5(k') ~ |(5(k)p, which is a real number. 
The second statistic of interest is the bispectrum B, defined by, 

(5(ki)<5(k2)<5(k3)) = (2^)35^(ki + k2 + k3)i?(ki, k2, kg). (4.2) 

The Dirac delta function , ensures that the bispectrum is defined only for k-vector configurations 
that form closed triangles: J2i = 0- Note that once the average is taken, the imaginary part of 
(5(ki)i5(k2)(5(k3) goes to zero. 

It is convenient to define the reduced bispectrum Q123 = (5(ki,k2,k3) as, 

O ^ B{k„k,,^,) 

"^''^ - P{k,)P{k2) + P{ki)P{ks) + P{k,)P{k3) ' ^ ■ ' 

which takes away most of the dependence on scale and cosmology. The reduced bispectrum is useful 
when comparing different models, since it has a weak dependence on cosmology and one can thus 
break degeneracies between cosmological parameters to isolate the effects of gravity. Hereafter, when 
we speak of the bispectrum we are always referring to the reduced bispectrum. 



5 Results 



In this paper, we present two ways of comparing the f{R) and ACDM reduced bispectra we obtain 
from N-body simulations. The differences depend on whether the models are matched in their initial 
or final power spectra. In method A we compare the output bispectra from N-body simulations with 
the same initial power spectra. Thus some of the difference in the bispectra can be attributed to 
the different amounts of final nonlinear power in the two sets. Method B tries to separate these 
contributions by generating modified initial power spectrum ACDM simulations whose power spectra 
at z = match those of the f{R) simulations. 

For both methods we compute the bispectrum randomly drawing fc-vectors from a specified bin, 
namely AA; and randomly orientating the triangle in space. We make the number of random triangles 
to depend on the number of fundamental triangle per bin, that scales as fcifc2^3Afc'^ [14]. In this 
paper we always choose Ak — 3k,nin- For the equilateral case, at scales of fc ~ 0.65 /i/Mpc we are 
generating ~ 5 x 10* triangles. We have verified that increasing the number of triangles beyond this 
value does not have any effect on the measurement. 

5.1 Method A (matched initial power spectrum) 

Our first test of f{R) vs ACDM bispectra utilize the same initial power spectrum. In our f{R) models, 
modifications to gravity go to zero rapidly with redshift and the expansion history differs negligibly 
from ACDM. Thus models with the same initial power spectra as ACDM fit observations at high 
redshift, such as primary CMB anisotropy, equally well. 

Since all N-body simulations start from the same initial power spectrum, the f{R) modifications 
to gravity during the acceleration epoch lead to differences in the dark matter power spectra at low 
redshift that increase with |/flo| as was noted in Fig. 2 of [10]: these differences reach up to ^ 50% 
for l/^ol = lO"** and - 10% for J/floj = 10"^ for k ~ 0.5/i/Mpc with respect to the ACDM model. 

Bispectra for matched initial power spectra but differing final power spectra is what is usually 
computed analytically [1, 2]: one predicts (using the modified Euler and continuity equations in 
perturbation theory) the reduced bispectra for different gravity models starting from a given initial Sk 
field. One might expect that the reduced bispectra differences are independent of the power spectra, 
but this is only strictly true for equilateral configuration and only in the tree-level regime (for k < 0.06 
h/Mpc at z = 0). This is the main caveat of method A: the differences seen in the reduced bispectrum 
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Figure 1. Relative dark matter reduced bispectrum deviations (following method A) between ACDM and 
f{R) models for k2 — 2ki = 0.4 h/Mpc (left panels) and equilateral configuration (right panels) at 2: = 
as a function of the angle between ki and k2, namely 812 (left panel) and as a function of k (right panel) 
for l/iiol = 10~*, 10~^, 10~® (top to bottom). Blue points (squares) correspond to chameleon simulations 
and red points (circles) to non-chameleon. Both ACDM and f{R) bispectra have been computed from N- 
body simulations with the same initial conditions. As a consequence, the corresponding final {z — 0) power 
spectra of the compared models are different. Error bars are the 1-a standard deviation of the ratio of Q 
values amongst the 6 independent runs. Because of that, the errors due to cosmic variance cancel out. Only 
L = 400 Mpc/h side-box runs are used. 

could be due to differences in the final matter power spectrum and not unique signatures of f{R) 
gravity. 

In Fig. 1 we show the dark matter reduced bispectra deviation between ACDM and f{R) for 
^2 = 2fci = 0.4/i/Mpc (left panels) and for equilateral configurations (right panels) for 2; = according 
to method A. The top panels correspond to f{R) theories with l/^ol = 10~^; |/flo| — 10~^ for middle 
panels; and |/flo| = 10^^ for bottom panels. The blue points correspond to full f{R) theories whereas 
the red points to non-chameleon ones. Deviations of f{R) bispectra with \fiio\ ~ 10^* with respect 
to ACDM present a characteristic shape dependence, where the difference is maximal for 6*12 ~ and 
TT, and minimal for 612 ^ O.Gtt, for both chameleon and non-chameleon and it increases as the scale 
is reduced. A similar trend is present for non-chameleon theories with |/ijo| = 10^^. On the other 
hand, chameleon theories with \fRo\ = 10~^ present a constant deviation from ACDM of ^ 10%. 
For l/jjol = 10~^, both chameleon and non-chameleon present a constant deviation from ACDM of 
< 5% and is consistent with 0. The errors of Fig. 1 are suppressed compared to the individual cosmic 
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variance errors because we are taking the ratio of N-body simulations with the same initial power 
spectrum and phases. Therefore, we can conclude that having the same initial conditions, the dark 
matter bispectra of ACDM and f{R) theories is significantly different for \frio\ ^ 10^^, especially for 
elongated triangles {612 — 0, tt) and can reach deviations in the reduced bispectra up to ~ 10% for 
l/flol = 10~^ and up to ^ 12% for I/aoI = 10~^. We see a similar dependence on triangle shape as 
shown in Fig. 5 of [1] (note that /3 as defined there is l/\/6 for /(i?)). 

Although differences in the final dark matter power spectra between ACDM and f{R) theories 
of the same initial power are large and potentially easier to test than those in the bispectra, it is 
possible that the galaxy power spectra for ACDM and f{R) models could still be similar for some 
particular galaxy bias model [18]. Since the galaxy bias acts differently on the power spectrum and 
on the bispectrum, it would be very unlikely that the same galaxy bias could make Pgt?°^ = Fg^f ^ 

and Q^^i^^ = Q^if simultaneously. 

Conversely, changes in the initial power spectra between the models might conspire to make an 
f{R) model look like a ACDM model for the power spectrum at z = 0. These can be hidden from 
the CMB at high redshift if they only occur at high k. Because of that, in the next section we assess 
the differences between the ACDM and f{R) dark matter reduced bispectra in the case where both 
models have the same final power spectrum. 

5.2 Method B (matched final power spectrum) 

The final power spectra of the f{R) models deviate significantly from that of the ACDM model with 
the same initial conditions (see Fig. 2 of [10]). These deviations reach ~ 50% for \fRo\ = 10^"*; ^ 10% 
for |/_Ro| = 10^^; all at fc ~ 1 h/Mpc and at z = 0. That begs the question of whether bispectrum 
differences seen in Method A are driven by these final power spectrum differences or by uniquely 
gravitational modifications. 

To address this question, we would like to adjust the initial conditions of the f{R) simulations 
until the final power spectra match that of ACDM at 2; = 0. However, the f{R) simulations are com- 
putationally very expensive (a factor of ^20 increase over ordinary GR simulations). Instead we do the 
converse: we adjust the initial conditions of the ACDM model until its final power spectrum matches 
the f(R) simulations. Matching ACDM to the f{R) simulations still tests whether the remaining bis- 
pectra difference between the two models reflects gravitational modifications, independently of power 
spectrum differences. 

5.2.1 Power Spectra Matching 

In order to match final power spectra we need a means of quickly predicting the impact of adjusting 
initial conditions in ACDM. HaloFit [17] provides an approximate analytic mapping between the 
initial and final power spectra. Using HaloFit we can determine the desired initial conditions, run the 
matching ACDM simulations and compare the bispectra with those of the f{R) simulations.^ 

We first test the accuracy of HaloFit in modeling the ACDM simulation results (see Fig. 2). 
In the HaloFit computation we use the same transfer function as employed in the ENZO code. We 
see that for the L = 400 Mpc/h runs the data-points with fc > A:jv/2 ~ 0.50 h/Mpc underestimate 
the value of the predicted power spectrum by HaloFit. Here, fc^r is the Nyquist mode defined as 

1 /3 

kjsi — nNp /(2L). We can in principle solve this limitation by using smaller boxes. For the L = 256 
Mpc/h runs, fcAr/2 ~ 0.79/i/Mpc and we see that up to this scale the simulation agrees with the 
theoretical prediction. However the errors increase considerably as we reduce the box-size. Error 
bars in Fig. 2 correspond to 1-a standard deviation amongst 6 independent runs. Likewise at low fc, 
the simulations carry large sampling errors even for the largest boxes. To evade these problems, we 

"^An alternative (much faster) approach would be, instead of re-running N-body simulations and computing the 
evolved bispectrum from there, to use fitting formulae from the literature to predict the ACDM bispectrum from the 
matched power spectrum which add a small running of the tilt. We have attempted this route; unfortunately we have 
found that available fitting formulae are not sufficiently accurate for our purposes [5] ; as we will see we need a relative 
accuracy of better than 4% on the reduced bispectrum which exceeds that of available fitting formulae in the mildly 
non-linear regime of interest. 
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Figure 2. HaloFit accuracy for the power law ACDM simulations with box size L = 256 Mpc/h (blue 
squares) and L = 400 Mpc//i box size (red circles) Left panel: linear power spectrum (dashed line) and 
non-linear prediction from HaloFit (dotted line) plotted with the simulation results. Right panel: fractional 
difference of simulation results and the HaloFit prediction. In all cases the errors correspond to 1-a standard 
deviation amongst 6 independent runs. 



use HaloFit to model only relative differences between simulations of the same L = 400 Mpc/h size, 
resolution and initial phases as we shall now describe. 

In order to match the excess small scale final power in the /(i?) model, we add an extra running 
of the spectral tilt parameter to the ACDM initial power spectrum. Specifically, we assume a 3 
free-parameter initial power spectrum model: 

^.(^)=A(^j , (5.1) 

where Pq is the amplitude of the power spectrum at kp ~ O.l/iMpc^^ and z = without the effect of 
the transfer function and a is the running of the tilt, 

=no + ah-i{k/kp). (5.2) 

We therefore have 3 parameters p — {Pq, uq, a} which specify the initial conditions. 

To find the best-fitting three parameters for a given model, we take the simulation results for 
the power spectrum ratio (following method A), 

«.Ut)^l0^\ (5.3) 

Next we use the HaloFit prescription P^^(fc; p™'^*'^^) for the non-linear matter power spectrum at 
z = to find the best parameter set p™^*'^^^ minimizing the given by 



fHaloFit(fcj;P°) 



(5.4) 



where the sum runs over bins in k. Here, p" describes the initial power spectrum used for the f{R) 
simulations (see first line in Tab. 1). Finally we simulate a matched ACDM simulation with the same 
initial phases as the original but with a rescaling of the initial power 

P (k - n™atch\ 

P,c{k;p-^'^^) = Piao.ig.(fc) '),^(P.pO) ^ (5.5) 

In order to avoid confusion, we designate these ACDM simulations as "matched"; whereas ACDM 
without this modifier denotes the standard, power law, initial conditions (the one used in §5.1). 
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Table 1. Best-fit linear ACDM initial power spectrum parameters Po, no, a (see text) tliat matcli tlie N-body 
power spectrum at z = for each f{R) simulation. Tliese values have been computed minimizing the of the 
N-body non-linear power spectrum and the non-linear HaloFit power spectrum for k < fc]v/2 ~ 0.5 /i/Mpc. 
The same transfer function and the same cosmology is used in all cases. Also ag is shown for clarity. 
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Figure 3. Relative power spectrum offset between f{R) simulations and corresponding matched-ACDM 
simulations for different |/ho| values: for the chameleon simulations with I/hoI = 10"' * (red circles), IQ-^ 
(blue triangles), 10"® (green stars); and for the non-chameleon simulations with j/jjo| = 10"^ (pink diamonds), 
IQ-^ (orange hexagons), IQ-^ (cyan pentagons). The matched simulations make use of initial power spectrum 
conditions shown in Table 1 found by fitting relative deviations with HaloFit. 



The advantage of this matching method is that we only model relative deviations with the HaloFit 
prescription. Thus the cosmic variance of the original simulations scale out as do absolute errors in 
HaloFit, initial condition generators, resolution, etc. In Table 1 we show the best-fit values of the 3 
initial power spectrum parameters. We have only used i?sim(^ ^ 0-5 /i/Mpc) for the minimization. 

In Fig. 3 we show Pf (r) / Pmatchcd — Ij where Pf(R) is the power spectrum of the f{R) simulations 
as before and Pmatchod is the power spectrum of the matched ACDM simulations. Fig. 3 indicates 
that HaloFit is an excellent tool to predict relative differences in non-linear power spectra even for 
(some) non-standard ACDM models. We see that the differences between the matched ACDM and 
f{R) power spectra are up to ^ 4% in the range 0.1/i/Mpc < k < lft,/Mpc, although for most of the 
scales and the cases are about 2-3%. 

As an aside, we can also test the absolute accuracy of HaloFit's prediction for the power spectra 
of the matched models. Examples for different matched models are shown in Fig. 4. HaloFit produces 
a good fit compared with the sample variance errors for all k < 0.5Mpc//i. As in the pure ACDM 
case, the sample variance at low k in the simulations is quite large. Deviations up to 10% for k < 1.0 
Mpc//i likewise appear due to the limited simulation resolution. Our modeling of relative effects 
eliminates these small differences. 
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Figure 4. HaloFit accuracy for a representative set of matciied ACDM models (matclied to chameleon 
/flo| ~ 10~*, black squares; 10~^ red circles; 10~^, blue triangles). HaloFit is accurate within the errors for 
all k < 0.5Mpc//i for these models which contain running of the tilt. Deviations from sample variance and 
resolution in the simulation seen here are largely absent in the relative matching technique shown in Fig. 3 



In reality one does not observe at a single z but in a wide z-range. As mentioned above it is not 
possible to match the power spectrum at widely separated redshifts simuhaneously and this feature 
can provide observational signatures independent from the bispectrum. We can quantify this further 
by estimating over what redshift interval the power spectrum matching is expected to hold. Changes 
in the Pf(R) / Pacdm were studied in detail by [10] and the excess evolves on the Hubble time scale. 
Therefore we generically expect that the matching evolves across a redshift interval of Az ~ 1, i.e. 
no faster than any other aspect of the modeling. 

5.2.2 Bispectrum 

With the simulations of the matched ACDM models, we can now compare the bispectra for ACDM 
and f{R) models whose final power spectra match to a few percent. 

In Fig. 5 we show Q/(i^)(fc)/Qmatchcd — 1 for fc2 = 2fci = 0.4/i/Mpc (left panel) and for equilateral 
triangle configuration (right panel), where Qf(R){k) is the reduced bispectrum for f{R) simulations, 
and Qmatchod is the reduced bispectrum for the matched ACDM simulations. Red points show the 
ratio for non-chameleon simulations whereas blue points for chameleons ones. Top panels correspond 
to |/_Ro| — 10""*, middle panels to \fm\ — and bottom panels to \fRo\ — 10^^. In particular, 

we see that for the chameleon and non-chameleon cases with |/flo| = 10~^ and 10^^ the deviation is 
very close to (< 2%). For the |/ko| = 10^^ some differences appear: for the non-chameleon case 
there is an excess of ~ 4% and for the chameleon case there is a deficit of ^ 4% in Qf(R) respect 
to Qmatchod, both within 5-6(7. The value |/flo| ^ 10~^ is special in that it marks the onset of the 
chameleon mechanism in the largest structures in the simulations. The chameleon effect may have 
a small but measurable impact on Q in this transition region where the chameleon effect is present 
for some but not all structures. Analogous transient enhancements appear in the mass function [7]. 
One should bear in mind though that this difference is of order the difference in the matched power 
spectra which varies between the full and no-chameleon cases. 

Thus, for all values of /j^q deviations are below ^ 4%. In particular we do not observe that 
squeezed triangles (those with 612 ~ 0, tt) present higher deviations between different gravity models 
as has been observed in method A (Fig. 1) and predicted from theoretical models that followed the 
same assumptions as adopted in method A [1, 2]. 

Finally, we found that it is better to analyze the deviation between reduced bispectra Q rather 
than between bispectra B. This is because the power spectrum dependence is partially canceled in 
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Figure 5. Relative reduced bispectrum deviation for matched final power spectra (method B) between 
f{R) and ACDM simulations for ~ 2fci — 0.4/i/Mpc as a function of 612 (left panels), and equilateral 
configuration as a function of k (right panels), all at z = 0. Upper panels correspond to |/ijo| = 10"**, 
middle panels to l/ijol = 10~^ and bottom panels to \fRo\ ~ 10^^. Red points correspond to non-chameleon 
simulations, whereas blue points to full f{R) simulations. The error-bars are the la standard deviation 
amongst the ratio of 6 independent runs. Since we are taking the ratio between runs with the same initial 
phases, the cosmic variance errors are not present. 



the reduced bispectra. In spite of having run ACDM simulations to match the f{R) power spectra, 
several percent differences between f{H) and matched ACDM power spectra are still present (Fig. 3). 
These lead to higher deviations in B between the models (up to ~ 8% in some cases) than in Q. 
Thus, using Q instead of B is much more robust if we want to compare models with similar power 
spectra. Of course, one should keep in mind that not all the P(fc)-dependence is cancelled when using 
Q, as evidenced by comparing with the results of method A. Strictly speaking, this is only true for 
equilateral configurations and up to tree-level in Eulerian perturbation theory. 

Finally, one may want to make a connection between these results and some analytic model, 
namely perturbation theory (PT). Since at tree level in PT the reduced bispectrum is independent of 
the power spectrum (at least for equilateral configuration), the differences observed between Fig. 1 
and 5 should be due to higher order corrections in ACDM. At 1-loop, corrections to the bispectrum 
can be found in e.g., [14, 16]. One can see that the leading terms depend on the linear and one-loop 
power spectrum and weakly on cosmology and gravity through the standard tree level bispectrum 
kernel (see [1] for a modification of this kernel for some f{R) theories). Thus, a small modification of 
this formula could be expected due to f{R) gravity. However, this interpretation should be considered 
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more in a qualitative way than in a strictly quantitative way. In fact one should take into account that 
the precision of 1-loop PT for the bispectrum is not much better than the other (phenomenological) 
analytic formulae [5] . As we have already mentioned, currently there is no analytic model that predicts 
the bispectrum at scales of interest here with an accuracy of few percent. 

5.2.3 Discussion 

If the remaining ~ 4% deviation for |/i?o| = 10~^ reflects gravity and not the residual mismatch 
in power spectra, then it is in principle measurable with large-volume surveys. In this work, con- 
sidering only the 6 runs of 400 Mpc/h box-size and provided that h — 0.73, the total volume 
is 6 X (OAGpc/h)^ ~ IGpc'^. We expect that future surveys will cover larger volumes: BOSS'^ 
V - 5(Gpc//i)3, DES'^ V - 10(Gpc/h)3 or EUCLID V ~ 100 (Gpc//i)3. As the 6 runs have dif- 
ferent initial conditions we can use them to estimate the expected error on Q in the limit that it is 
dominated by cosmic variance. We have measured that the error in Q for our simulations at scales of 
k ^ 0.3/i/Mpc is about 5%. We assume that the variance scales as the inverse of the number of modes, 
and thus the standard deviation approximately scales as V^^^ . Therefore, for a 10 Gpc"^ survey the 
error bars could, in principle, be as much as \/lO ^ 3 times smaller than our prediction. This implies 
that a survey with > lOGpc^ volume (e.g., DES, EUCLID) would yield an error on Q ^ 2% at these 
scales. Since the expected deviation may be of order 4%, having smaller errors would help us to 
confirm or discard possible deviation of the bispectrum due to modifications of gravity. 

On the other hand, we have analyzed the dark matter bispectrum which is not directly observable. 
In practice, sources of error that we have neglected here may appear: i) galaxy-surveys provide a biased 
information about dark matter, ii) additional effects such as redshift distortions change the observed 
bispectrum (in fact we expect modified gravity to affect redshift distortions more strongly than the 
density field itself). Also as we go to higher z, we expect less deviations at a given scale. Conversely, 
the matched power spectra at z = would become mismatched and provide other observable effects. 

The results from Fig. 5 provide another important result. We have seen that two f{R) theories of 
gravity with indistinguishable non-linear dark matter power spectrum, have very similar and possibly 
indistinguishable dark matter bispectra. This opens up the possibility of using these two statistics 
to break degeneracies in the galaxy bias in a way that is robust to the assumptions about the true 
underlying gravity model. 

In fact the f{R) effects on the power spectrum are at the 20-50% level. A modification of galaxy 
bias achieving similar effects would likely affect the reduced bispectrum at least at the 10% level (for 
example, a linear bias term affect the power spectrum cx 6^ and the reduced bispectrum cx l/6i), 
significantly larger than the f{R) effects on the reduced bispectrum. 

6 Conclusions 

In this work we have analyzed the deviations in the reduced bispectrum produced by a modification 
of gravity, specifically the f{R) class of models, both with and without the chameleon mechanism. In 
order to do that, we make use of a suite of f{R) and ACDM simulations. We have proceeded in two 
different ways to analyze the bispectrum deviation from these simulations, methods A and B, which 
differ in whether the initial or final power spectra of the two cosmologies are set equal. 

Method A compares the bispectrum output of f{R) and ACDM N-body simulations with the 
same initial power spectrum. Fig. 1 shows the bispectrum deviation obtained using this method. 
We observe a considerable deviation (up to 10 — 15%) in the reduced bispectrum between these f{R) 
models and the ACDM one. Such differences in the bispectrum could be easily detected by surveys 
covering volumes > 1 Gpc such as e.g., the on-going BOSS survey. Higher deviations are seen for 
higher values of Ifm] and for squeezed triangle configurations. In this method, both ACDM and 
f{R) gravity runs start from the same initial values. Because of that, the different evolution of 
the gravity models naturally leads to different power spectra (as was observed in [10]) and also to 
different bispectra. This way of proceeding is equivalent to the theoretical works of Bernardeau & 

^Baryon Oscillation Spectroscopic Survey 
^Dark Energy Survey 
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Brax [1] , Borisov & Jain [2] . In order to explain discrepancies between the matter power spectrum in 
f{R) and the observed galaxy power spectrmxi, one could invoke a scale-dependent galaxy bias. Since 
galaxy bias enters into the reduced galaxy bispectrum in a different way than in the power spectrum, 
bispectrum measurements can in principle close this loophole. 

Alternatively, the large power spectrum differences can be eliminated by changing the shape of 
the initial power spectrum to instead match the final dark matter power spectrum at z = 0. This is 
at the base of method B. In this method, we compute the bispectrum deviation between a ACDM 
and a f{R) model, both with the same final power spectra. Thus, we simulate a ACDM model 
with certain initial power spectrum parameters (summarized in Table 1) that are adjusted to best 
match the /(i?) power spectrum at z = 0. From the simulations outputs we compute power spectra 
and bispectra. For the power spectra, residual differences are never higher than 4% in the range 
O.l/i/Mpc < fc < 1/i/Mpc. 

Likewise the differences in the reduced bispectrum are also smaller in the matched comparison. 
For the |/flo| = 10""* and 10"^ cases, the Q deviation is consistent with within la. For |/flo| — 10 
deviations in Q at most reach the 4% level with 5 — 6a significance. These deviations are potentially a 
signature of the onset of the chameleon mechanism in the largest structures in the Universe. However 
given that this is the same order as the power spectrum difference it is unclear whether these differences 
indicate power-spectrum-independent modified gravity effects or that the two power spectra are not 
perfectly matched. In the former case, larger surveys like EUCLID will allow for a measurement of 
the bispectrum with enough precision to obtain a > 6cr significance, even when exactly matching the 
power spectra. 

On the other hand, the effect of deviations from GR gravity on the reduced bispectrum are weak 
compared to those on the power spectrum (at least for the cases considered here), opening up the 
possibility of breaking the galaxy-bias degeneracy. In fact the effect of galaxy bias is expected to be 
different in the power spectrum and in the bispectrum, which is why, in the context of GR gravity, 
the bispectrum is used to constrain galaxy bias. While the shape of the non-linear power spectrum 
seems to carry information about the underlying gravity model, one may always argue that a shape 
of the evolved power spectrum not compatible with the GR predictions could be due to biasing. For 
the cases we have considered here, the dependence of the reduced bispectrum on deviations from 
GR is weaker than the effects of bias modifications necessary to explain the deviations in the power 
spectrum. While we have only studied f{R) models here, there is no apparent reason why this result 
should be specific to f{R)- Hence, if our findings were to remain qualitatively true for other gravity 
modifications, this would confirm the usefulness of employing the reduced bispectrum together with 
the power spectrum to constrain bias parameters. 
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